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ABSTRACT 

We present a numerical method for spatially 1.5-dimensional and time- 
dependent studies of accretion disks around black holes, that is originated from 
a combination of the standard pseudo-spectral method and the adaptive domain 
decomposition method existing in the literature, but with a number of improve- 
ments in both the numerical and physical senses. In particular, we introduce a 
new treatment for the connection at the interfaces of decomposed subdomains, 
construct an adaptive function for the mapping between the Chebyshev-Gauss- 
Lobatto collocation points and the physical collocation points in each subdomain, 
and modify the over-simplified 1-dimensional basic equations of accretion flows to 
account for the effects of viscous stresses in both the azimuthal and radial direc- 
tions. Our method is verified by reproducing the best results obtained previously 
by Szuszkiewicz & Miller on the limit-cycle behavior of thermally unstable ac- 
cretion disks with moderate viscosity. A new finding is that, according to our 
computations, the Bernoulli function of the matter in such disks is always and 
everywhere negative, so that outflows are unlikely to originate from these disks. 
We are encouraged to study the more difficult case of thermally unstable accre- 
tion disks with strong viscosity, and wish to report our results in a subsequent 
paper. 

Subject headings: accretion, accretion disks — black hole physics — hydrody- 
namics — instabilities 
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Introduction 



The radiation pressure-supported inner region of geom etrically thin, optically th ick 
Shakura-Sunyaev accretion disks (SS D) around black h oles (IShakura &: SunyaevI Il973l ) is 
known to be thermally unstable (e.g. iKato et al.l Il998l . Chap. 4), but the occurrence of 
an instability does not necessarily mean that the disk will be disrupted after the charac- 
teristic growth-time. A possible fate of the thermally unstable inner region of SSDs is the 
so-called limit-cycle behavior, i.e., the nonlinear oscillation between two stable states. Simi- 
lar to the cas e of dwarf novae, the limit-cycle behavior was realized from a local and steady 
analysis (e.g. iKato et al.lll998l . Chap. 5), i.e., from an S-shaped sequence of steady state so- 
lutions at a certain radius in the M — S (mass accretion rate vs. surface density) plane, 
with the lower and middle branches of the S-shaped sequence corresponding to stable gas 
pressure-supported SSD solutions and unstable radiation pressure-supported SSD solutions, 
respectively, and the upp er branch corresponding to stable slim disk solutions constructed by 
Abramowicz et al.l (Il988l ): and has been jus tified by a number o f works performing global and 



time-dependent numerical con ip utations jHonma et al. 1991 : Szuszkiewicz fc Miller 1997 



19981 . 12OO1I : IXeresi et al.l l2004al lbl: iMaver fc Prindel bood ). Unlike the case of dwarf novae 
however, only one astrophysical object, the Galactic microquasar GRS 191 5-1-105, has been 



known to show the theor e tically predicted limit-cyclic luminosity va. r iations (INavakshin et al. 



200a IJaniuk et al.ll2002l : IWatarai fc Mineshigell2003l : IOhsugall2006l : iKawata et al.ll2006h . 



We select the paper of ISzuszkiewicz &: Millerl (120011 . hereafter SMOl) as the representa- 
tive of existing theoretical works on the limit-cycle behavior of black hole accretion disks with 
the following two reasons. First, SMOl adopted a diffusion- type prescription for viscosity, 
i.e., the rcf) component of the viscous stress tensor is expressed as 

dVL 



Tr4, = aHcspr 



dr ' 



1) 



where p is the density, Q is the angular velocity, Cg is the sound speed, H is the half-thickness 
of the disk, and a is a dimensionless constant parameter; whereas all other relevant works 
used a simple prescription 

rr<j, = -ap, (2) 

where p is the pressure, and a is also a dimensionless constant but has been rescaled (denoting 
a in expressions [1] and [2] as ai and 02, respectively, then a2 = [3\/6/2]ai). It is known 
that the direct integration of the differential equations describing transonic accretion disks 
with the diffusive form of viscosity is extremely difficult, while that with the ap viscosity 
prescription becomes much easier (see the discussion in SMOl). It should be noted, however, 
that expression ([2]) is only an approximation of expression ([1]) under a number of conditions 
(including assuming that the disk is stationary, geometrically thin, Newtonian Keplerian 
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rotating, and in vertical hydr ostatic equilibrium, e.g . Kato et al.l 119981 . Chap. 3). More 
seriously, as shown recently by lBecker fc SubramanianI (120051 ). expression ([1]) is the only one 
proposed so far that is physically consistent close to the black hole event horizon because 
of its diffusive nature, whereas expression ([2]) as well as some other viscosity prescriptions 
would imply an unphysical structure in the inner region of black hole accretion disks. Second, 
SMOl did complete very nice numerical computations, all the curves in their figures showing 
the evolution of disk structure are perfectly continuous and well-resolved on the grid; while 
some fluctuations appear on the curves in the figures of other relevant works, which might 
make one to worry whether there had been some hidden numerical instabilities in the code. 

As evidenced by SMOl, thermally unstable accretion disks undergo limit-cycles when 
viscosity is moderate, i.e., the viscosity parameter a ~ 0.1 (hereafter all the numerical values 
of a are for a2 unless otherwise specified); and the instability seems to be catastrophic when 
viscosity is weak, i.e., a ~ .001. On the other hand, in the case of very strong viscosity, 
i.e., a ~ 1. IChen et al.l (119951 ) found that the S-shaped sequence of steady state solutions in 
the M — S plane does not form, instead, sli m disk solutioris and optically thin advection- 



dominated accretion flow (ADAF) solutions (iNarayan fc Yi 



1994 



Abramowicz et al.lll995r ) 



are combined into a single straight line. Accordingly, Takeuchi fc Mineshigel ( 1998 ) per 



formed time-evolutionary computations using the viscosity prescription with a = 1 and 
proposed another possible fate of thermally unstable accretion disks: the very inner region of 
the disk finally becomes to be an ADAF-like feature, while the outer region keeps being the 
SSD state, forming a persistent two-phased structure. While this result is really interesting 
since a phenomenological SSD-I-ADAF model has been quite suc cessfully applied to black 
hole X-ray binaries and galactic nuclei (e.g., iNarayan et al.lll998l ). SMOl stated that they 
could not make computations for a = 1 because of difficulties in keeping their code numeri- 
cally stable, and pointed out that i t is wo rth checking whether the persistent ADAF feature 
obtained in iTakeuchi fc Mineshigd (119981 ) would survive changing the viscosity prescription 
to the diffusive form. 

We purpose to study thermally unstable accretion disks if they are not disrupted by 
instabilities, that is, we wish to check whether the limit-cycle behavior is the only possible 
fate of these disks provided viscosity is not too weak, or a transition from the SSD state to 
the ADAF state is the alternative. As in SMOl, we adopt the diffusive viscosity prescription 
of equation ([1]) and make spatially 1.5-dimensional, time-dependent computations. But we 
choose a numerical method that is different from either of SMOl or of iTakeuchi fc Mineshige 
(119981 ) ■ and that is the adaptive pseudo-spectral domain decomposition method. With this 
method, we hope to be able to perform computations for various values of a ranging from 
~ 0.1 to ~ 1, and to obtain numerical results at the quality level of SMOl. In this paper, 
we describe our numerical algorithm and techniques in details and present computational 
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results for a = 0.1 as a test of our algorithm. We wish to report our resuhs for larger values 
of a in a subsequent paper. 



2. Numerical Algorithm 

As the main intention of this paper, in this section we present a numerical algorithm to 
solve a partial differential equation (or equations) in the general form 

— — = L{u{r,t)), r e [rmin,rmax], (3) 

where u{r,t) is a physical quantity that is a function of the spatial independent variable 
r (e.g., the radius in the cylindrical coordinate system) and the time t, and L is a partial 
differential operator of r and can be linear or nonlinear. 



2.1. Scheme of Spacial Discretization 



We first describe the standard Chebyshev pseudo-spectral method that is used to dis- 
cretize the spatial differential operator L. This method has been explained in several text- 
books ( 



Gott 



ieb fc Orszad Il983l : ICanuto et all Il988l : iBovdl l2000l : iPevretl 120021 ) . Recently, 



Chan et al.l (l2005l . 120061 ) applied it to studies of astrophysical accretion flows and discussed 



its advantages. 

Concretely, a series with finite terms is used to approximate a physical quantity u{r) as 



uirk) = u[g{rk)] 



N N /nk7r\ 

UnTn{fk) =^Un COS i -j^ j 
n=0 n=0 ^ ^ 



(4) 



where Tn{fk) is the n-th order Chebyshev polynomial; [k = 0, 1, 2, ...N) is the Chebyshev- 
Gauss-Lobatto collocation points and is defined as fk = cos{kn/N), with N being the number 
of collocation points; = g{fk) is the mapping from the Chebyshev- Gauss-Lobatto colloca- 
tion points fk G [—1, 1] to the physical collocation points G [rmin,rmax] that is a strictly 
increasing function and satisfies both g{—l) = and g{l) = r^ax] Un is the spectral 
coefficients and can be calcu lated from the ph ysical values u{rk) by a fast discrete cosine 
transform (hereafter FDCT, iPress et al.l Il992l . Chap. 12); contrarily, if one has then 
u{rk) can be obtained immediately by a inverted FDCT. 



The radial derivative du{r)/dr is also a function of r and in principle can also be 
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approximated by a series that is obtained by using the chain rule 

du{rk) 1 du[g{fk)] 1 ^ 



^u^Tn{rk). (5) 



dr dg/dr dr dg/dr 

The spectral coefficients can be calculated from by a three-term recursive relation 

u'n = 0, 
u'n-i = '2Nun, 

Cnu'n = Mi,+2 + 2(n + l)n„+i, (6) 

where Cq = 2, and c„ = 1 for = 1,2, A^. Subsequently, 9M[5'(ffc)]/9f is calculated from 
u'^ by a inverted FDCT, and then substituted into equation ([5]) to obtain discrete spatial 
derivatives du{rk)/ dr. 

To summarize, we define a discretized differential operator D for the continuous differ- 
ential operator 8/ dr. The operator D carries out the following works: (1) using FDCT to 
calculate from u{rk); (2) using the three-term recursive relation equation to obtain 
from iin'-i (3) using a inverted FDCT and equation ([5]) to obtain du{rk)/dr. Finally, 
we use D to construct a discretized operator L to approximate the operator L in equation 
([3]). For example, if L = dr{uidrU2), where dr denotes d/dr, then L can be constructed as 
L = D[u,D{u2)]. 



2.2. Scheme of Time-Discretization 



We adopt two schemes to perform the time-integrati on, that is, we use a third order 
total variation diminishing (TVD) Runge-Kutta scheme (IShu &: Osherl Il988l ) to integrate 
the ffist two time-steps, and then change to a low C PU-consump tion scheme, the so-called 
third order backward-differentiation explicit scheme (jPeyretl |2002| . pp. 130-133), to carry out 
the rest computations. 

The third order TVD Runge-Kutta scheme is expressed as 

'(1) = M" + AtL(M"), 



n+l 



-U 
4 



4 4 ^ ^ 



u 



-u" + -u 



(2) 



-AtL(«(^)), (7) 

> of the physical quantity u at the 
and (n + l)-th time-levels, respectively; and u^^^ and m^^-* are two temporary variables. 



3 3 3 

where At is the time-step; and u""*"^ are the values of the physical quantity u at the n-th 
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The third order backward-differentiation exphcit scheme can be written as 

3 2 



where 



and 



3=0 3=0 



1 1 

ao = 1 + r- + 



ai 



a2 



1 -\- kn 1 + fcn + kn—1 
(l + /Cn)(l + fcn + fcn-l) 
^n(^n + ^n-l) 
1 + kji + kn—1 



knkn-li^ + kn 



^ ^ + kn 

kn-l{kn + + fcn + ^n-l) 



&0 



[l + kn){l + kn + kn-l] 

kn{kn ~\~ kfi-l) 
1 + /c„ + 



n^'n— 1 



kn — 



n— 1 



At ' 



with t", and being the times of the n-th, (n — l)-th, and {n — 2)-th time-levels, 
respectively. 

Of these two time-integration schemes, the former spends three times of the latter's 
CPU-time per time-step, but the latter is not able to start the time-integration by itself 
while the former is able to do. Therefore, we combine these two schemes in order to achieve 
a sufficient high order accuracy with minimal CPU-time consumption. 

Hereto, we have fully discretized equation In order to obtain a physically sound 
and numerically stable solution in a finite domain, it is additionally necessary to impose 
appropriate boundary conditions and to apply some filtering techniques to overcome the 
inevitable spurious nonlinear numerical instabilities in the code. We leave the details of 
these in the Appendix. 
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2.3. Domain Decomposition 



The numerical algorithm described in the above two subsections and the Appendix has 
been a usef ul implemeri t for so lving; partial differential equations and is essentially what was 
adopted in I Chan et al.l (120051 ) . However, it turns out that, as we have experienced in our 
computations, the above algorithm is insufficient for resolving the so-called stiff problem. 
This problem is a one whose solution is characterized by two or more space-scales and/or 
time-scales of different orders of magnitude. In the spatial case, common stiff problems in 
fluid mechanics are boundary layer, shear layer, viscous shock, interface, flame front, etc. In 
all these problems there exists a region (or exist regions) of small extent (with res pect to 



the g lobal characteristic length) in which the solution exhibits a very large variation (iPeyret 



f 



2OO2I . p. 298). When the Chebyshev pseudo-spectral method described in §2.11 is applied 



to a stiff problem, the accuracy of the method can be significantly degraded and there may 
appear spurious oscillations which can lead to nonlinear numerical i nstabilities or spuriou s 



predictions of solution behavior (the so-called Gibbs phenomenon, iGottlieb &: Shu 



1997|). 



The spectral filtering technique described in the Appendix is not able to completely remove 
these spurious oscillations, so that the solution is still not well resolved, and sometimes 
the computation can even be destroyed by the growing spurious oscillations. A special 
method that has been developed to overcome these difficul ties is the dornain de composition 
(IBayliss et al.l Il995l : iPeyretl |2002| ) . Here we mainly follow iBayliss et al.l (1X9951 ) to use this 
method, but with a different technique to connect the decomposed subdomains. 

The basic idea of domain decomposition is to divide a wide computational domain into 
a set of subdomains, so that each subdomain contains at most only one single region of rapid 
variation (i.e., with a stiff problem), and more grid points are collocated into this region by a 
special mapping function to enhance the resolution while the total consumption of CPU-time 
is not substantially increased. 

In each subdomain the solution is obtained by taking into account some connection 
conditions at the interface between the two conjoint subdomains. In general, appropriate 
connection conditions a re the continuities of both the solution and its spatial derivative 
normal to the interface (IBayliss et al.lll995l : iPeyretl |2002| ) . The continuity of the solution is 
satisfied naturally, but the continuity of its derivative cannot be achieved directly with the 
pseudo-spectral method because of the use of FDCT. To see this, let us divide the entire 



computational domain r G [r 



mini ' max_ 



into M subdomains. 



I L mini maxli 



where r, 



(1) 



^(1) 

' max 



mil 



4i) 



1,2,...,M, 



and ri^^^ 



(12) 

are the locations of 



(2) (2) _ (3) 
mm ■ rnm ; ' max — ^mm ' ^maa; — "^min ' ■ ' " '^'-^'^ ' max — ' max 

the interfaces between the subdomains. Because FDCT is used to calculate the numerical 
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derivative in each subdomain Si, one obtains two values of the derivative at each interface. 
Let d~u and d^u denote the left and right numerical derivatives of the physical quantity u 
at a certain interface, respectively, a seemingly rational choice for keeping the continuity of 
derivative is to set the numerical derivative at the interface to be the mean of d~u and d^u, 



I.e., 



(I) 

\ / interface 



d U + d^u 



(13) 



Unfortunately, in practice the connection technique of equation (fT3l) will often cause a nu- 
merical instability at the interfaces. 

We find that the connection between two certain subdomains Si and S'j+i can be nu- 
merically stable when their discretizations satisfy an additional practical condition. Let 



int\ 



' max 



denotes the location of interface between Si and Sj+i; rfl_^ and r\ 



N- 



1 ^ Si, r[ ^ ^ G Sj+i and r''^'_i < Tint < r\ 



denote the locations of the two nearest 



points to the interface, respectively; our computations show that if the condition 



r 



(i) 

N-1 



r 



1 



(14) 



is satisfied, then the connection of derivative represented by equation ( fT3i) will be numerically 
stable. 

If the stiff problem always appeared in a fixed spatial region, then the domain de- 
composition would be kept unchanged. However, in general this is not the c ase. Instead, 



the l ocation of region in which the stiff problem appears changes with time (iBayliss et al. 



19951 ). Therefore, the domain decomposition must be adjusted adaptively. To ensure the 
connection condition equation f|T^ at the interfaces of newly divided subdomains, an ad- 
justable mapping between the physical collocation points in each new subdomain Si and 
the Chebyshev-Gauss-Lobatto collocation points is needed. We adopt such a mapping in 
the form (see eq.[l]) 



9[r) 



— I r 



arctan 



vr 

a tan — r 
4 



r G [—1, 1] and r 



in the subdomain Si, which is a combination of two mapping functions. 



' max 



and 



arctan 



vr 



IT 



a tan — (r — 1) 
4 



Equation f|T6l) is a trivial linear mapp ing (IChan et al.ll2005l ). and equation 
ping proposed by iBayliss et al.l (119951 ). The parameter a in equation ([T 



(15) 
(16) 



(17) 

f|T71) is the map- 
is an adjustable 
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parameter, but equation f[T7j) is only a mapping from f G [—1, 1] to f G [—1, 1]. Therefore, 
we add equation f|T6|) in order to make a complete mapping from f G [—1, 1] to r^*) G Sj. 
The combined mapping equation f[T^ will concentrate the discrete grid points toward r^^-^ 
when a > 1 and toward rmaa; when a < 1, and will be reduced to equation ( fT6l) when a = 1. 



The adjustability of mapping equation (fT5|) is crucially important for achieving a nu- 
merically stable connection at the interfaces of subdomains. By substituting equation (fTSl) 
into equation (ITll) . we obtain 

(i+i) _ arctan [a^'^ tan f (n - 1)] } 



tan j{rr 



1) 



with 



U! = 



ri 



' max I min 
I max I rn.in 



COS , ^ 



rr 



cos 



(A^- l)7r 
N 



where a^*) and a^^^^^ are the mapping parameters for subdomains Si and Sj+i, respectively. 
Equation ( ITSl) can be used to determine the mapping parameters of every subdomain after 
giving a decomposition of computational domain {Si} and the mapping parameter a^^-* of 

]). As a result, we obtain a particular collocation 



mm: I max 



the innermost subdomain Si ( 
of discrete grid-points within the whole computational domain [rmin.rmax]- This colloca- 
tion ensures a stable connection of the derivatives between any two conjoint subdomains 
(eq.[T3]), and thus ensures a correct implementation of the pseudo-spectral method in each 
subdomain. The combination of the standard pseudo-spectral method and the adaptive do- 
main decomposition method finally names our numerical algorithm as that in the title of 
this paper. 



3. Limit-cycle Solutions 

We now verify our numerical algorithm by applying it to studies of thermally unstable 
black hole accretion disks with moderate viscosity and comparing our results with that of 
the representative work SMOl. 
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3.1. Basic Equations 

We write basic equations for viscous accretion flows around black lioles in the Eulerian 
form rather than Lagrangian form as in SMOl because partial differential equations in the 
Eulerian description take the general form of equation ([3]), to which our numerical algorithm 
is suited. The basic equations to be solved are 

_ _ dT, _T, d 

dt ^ dr r dr ^ ^ 

dvy ^ _^ ^ _ 1^ ^ - ^20) 
dt ^ dr p dr ' 

dl dl a d f ^ ^^^dil\ 



-VrTT + — r-^CsHE— , (21) 
dt or rS or \ or J 

dH dH 

dt ^ dr S (r — rg^ \ r J ^ 

dT _ dT 
dt ^ dr 



, T r aEcsHjrdQ/drf - p- _ 
12- 10.5/3 \ 0.67pif ^ 



H r dr 



(24) 



Equations flTI?]) . fl2U]) . fl^Tl) . and fl2^ are conservations of mass, radial momentum, an- 
gular momentum, and energy, respectively. As in SMOl, we adopt the diffusive form of 
viscosity (eq. [1]) in equations (pT!) and fl24|) : and abandon vertical hydrostatic equilibrium 
assumed in 1-dimensional studies, and instead introduce two new dynamical equations for 
the vertical acceleration (eq. [23]) and the evolution of the disk's thickness (eq. [22]), thus 
our studies can be regarded as 1.5-dimensional. In these equations f^, /, Ik, V^, Mbh, ^g, 
T, j3, and F~ are the radial velocity, speciflc angular momentum, Keplerian speciflc angular 
momentum, vertical velocity at the surface of the disk, black hole mass, gravitational radius 
(= 2GMbh/c'^), temperature, ratio of gas pressure to total pressure, and radiative flux per 
unit area away from the disk in the vertical direction, respectively. We use the 'one-zone' 
approximation of the vertically-averaged disk as in SMOl, so that Vr, fi, I, Ik, P, P, Cs, and 
T are all equatorial plane quantities, while Vz and F~ are quantities at the disk's surface. 
Additional deflnitions and relations of these quantities are 

P = Jj^ (25) 
/ I gMbht^ .... 



p 

= 1, (28) 

F- = ^±11 (29) 

3tr/2 + V3 + 1/tp' 

P = kpT + Prad, (30) 

F- f 2 \ 

Prad = — + ^ , (31) 



12c V VS, 

Tr = 0.34S(l + 6 X lO^V^-^'^), (32) 
P - " (34) 

P 

where equation ( l29l) is a bridging formula that is vahd for both optically thick and thin 
regimes, prad is the radiation pressure, and tr and Tp are the Rosseland and Planck mean 
optical depths. 



3.2. Specific Techniques 

As in SMOl, in our code the inner edge of the grid is set at r ~ 2.5rg, close enough to 
the central black hole so that the transonic point can be included in the solution; and the 
outer boundary is set at r ~ lO^rp, far enough away so that no perturbation from the inner 
regions could reach. A stationary transonic disk solution calculated with the ap viscosity 
prescription is used as the initial condition for the evolutionary computations. The ap 
viscosity prescription may seem inconsistent with the evolutionary computations adopting 
the diffusive viscosity prescription, but this does not matter. In fact, the initial condition 
affects only the first cycle of the obtained limit-cycle solutions, and all following cycles are 
nicely regular and repetitive. The time-step At is adjusted also in the same way as in SMOl 
to maintain numerical stability. We emphasize some techniques specific to our numerical 
algorithm below. 

The solution to be obtained covers so wide a range of the whole computational domain, 
and in particular, the thermal instability causes a violent variation of the solution (the 
stiff problem) in the inner region (inside 200rg) of the domain. In such circumstances the 
standard one-domain spectral method is certainly insufficient. Accordingly, we divide the 
whole computational domain into 6 subdomains and let each subdomain contain 65 grid 
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points, so the total number of grid points is 65 x 6 — 5 = 385 (there are 5 overlapping 
point at the interfaces of subdomains). At each time-level we apply the one-domain spectral 
method described in §2.11 into each subdomain. In doing this, various techniques are used 
to remove or restrain spurious non-linear oscillations and to treat properly the boundary 
conditions, in order to have a numerically stable solution in each subdomain as well as a 
stable connection at each interface. Then we use the scheme described in §2.21 to carry out 
the time-integration over the time-step At to reach the next time-level. 

In general, spurious oscillations are caused by three factors: the aliasing error, the Gibbs 
phenomenon, and the absence of viscous stress tensor components in the basic equations. The 
aliasing error is a numerical error specific to the spectral method when it is used to solve 
differential equations that contain non-linear terms, and c an be resolved by the spectral 



filtering techniqu e described in Appendix (see IPeyretl l2002l for a detailed explanation and 



Chan et al.ll2005l for a quite successful application of this technique). 



However, the spectral filtering itself cannot resolve the Gibbs phenomenon characteris- 
tic to the stiff problem, and the adaptive domain decomposition method described in §2.31 
becomes crucially important. In our computations, we set the spatial location where a large 
variation appears (e.g., the peak of the surface density S) as the interface of two certain 
subdomains and use the mapping equation f|T5l) along with the connection conditions f|T3l) 
and f|T^ . so that more grid points are concentrated on the two sides of the interface to 
enhance the resolution, and a stable connection between the two subdomains is realized. As 
the location of large variation is not fixed and instead shifts during time-evolution, we follow 
this location, redivide the computational domain and recalculate the mapping parameter 
for each new subdomain with equation (fTSj) . In practice, the stiff problem appears and its 
location shifts during the first about 30 seconds of each cycle whose whole length is ~ 700 
seconds. In these 30 seconds we have to redivide the domain and reset the grid after every 
interval of 0.1 seconds (or every 0.01 seconds for the first a few seconds), and the typical 
length of time-step At is about 10~^ — 10~^ seconds. For the rest time of a cycle (more than 
600 seconds), the stiff problem ceases, then the grid reset is not needed and the time-step 
can be much longer. 

In addition to the above two factors in the numerical sense, there is a factor in the 
physical sense that can also cause spurious oscillations. The viscous stress tensor has nine 
spatial components, but in 1-dimensional studies usually only the r0 component is included 
in the basic e quations, and omitting other components can result in numerical instabilities. 



In particular, ISzuszkiewicz fc Miller! (119971 ) noted that if the tensor components r^r and r^^ 
were neglected from the radial momentum equation, some instability would develop and 
cause termination of the computation because of the non-diffusive nature of the equation; 
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and that the instabihty was suppressed when a low level of numerical diffusion was added 
artificially into the equation. In our code we resolve the similar problem in a slightly different 
way. We add into the radial mo mentum equation (I2U1) t wo vi scous forces F^r and F^^ in the 



form as given by equation (5) of ISzuszkiewicz &: Millerl (119971 ). and accordingly add into the 



energy equation (12^ ) a heating term due t o visco us friction in the radial direction as given 
by equation (15) of ISzuszkiewicz fc Millen (119971 ). We find by numerical experiments that 



when a very small viscosity in the radial direction is introduced, i.e., with the radial viscosity 
parameter ~ 0.05a, where a is the viscosity parameter in the ap viscosity prescription, 
the spurious oscillations due to the absence of viscous stress tensor components disappear 
and the solution keeps nicely stable. 

Solving partial differential equations in a finite domain usually requires the Dirichlet 
boundary condition, i.e., changing the values of physical quantities at the boundary points 
to be their boundary values supplied a priori, or/and the Neumann boundary condition, i.e., 
adjusting the values of derivatives of physical quantities at the bounda r y poi nts to be their 



boundary values supplied a priori. For spectral methods, I Chan et al.l (120051 ) introduced a 
new treatment, namely the spatial filter as described in the Appendix, in order to ensure 
that the Dirichlet or/and the Neumann boundary condition can be imposed and numerical 
instabilities due to boundary conditions are avoided. We note, however, that their spatial 
filter treatment is applicable only for physical quantities whose boundary derivatives are 
equal, or approximately equal, to zero (e.g., the radial velocity Vr, its spatial derivative at 
the outer boundary is nearly zero); while for quantities whose boundary derivatives are not 
negligible (e.g., the specific angular momentum /, its spatial derivative at the outer boundary 
is not very small), the boundary treatment in finite difference methods still works, i.e., to 
let directly the values of those quantities at the boundary points be the specified boundary 
values. 

SMOl supplied boundary conditions at both the inner and outer edges of the grid, i.e., 
at the inner edge Dl/Dt = {d/dr + Vrd/dr)l = and dp/dr = 0, and at the outer edge 
both Vr and / are constant in time and / is slightly smaller than the corresponding Ik- In 
our computations we find that it is not necessary to supply a priori the two inner boundary 
conditions as in SMOl. With the two outer boundary conditions as in SMOl and free inner 
boundary conditions instead, we are able to obtain numerically stable solutions of the basic 
equations, and these solutions automatically lead to an almost zero viscous torque, i.e., 
Dl/Dt ~ 0, and a nearly vanishing pressure gradient, i.e., dp/dr ~ 0, in the innermost zone. 
This proves that the inner boundary conditions assessed in SMOl are correct, but we think 
that our practice is probably more natural and more physical. Once the state of an accretion 
flow at the outer boundary is set, the structure and evolution of the flow will be controlled 
by the background gravitational field of the central black hole, and the flow should adjust 
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itself in order to be accreted successfully into the black hole. In particular, both the viscous 
torque and the pressure gradient of the flow must vanish on the black hole horizon, and in 
the innermost zone a correct solution should asymptotically approach to such a state, thus 
in the computations no inner boundary conditions are repeatedly needed. 



3.3. Numerical Results 

We have performed computations for a model accretion disk with black hole mass 
Mbh = IOMq, initial accretion rate m = M/Mcr = 0.06 (M = — 27rrSfr is the accre- 
tion rate and M^r is the critical accretion rate corresponding to the Eddington luminosity), 
and viscosity parameter a = 0.1 (the equivalent a in the diffusive viscosity prescription is 
0.1 X (2/3^6) ~ 0.0272). It is known that the inner region of a stationary accretion disk 
with such physical parameters is radiation pressure-supported (/? < 0.4) and is therrn ally un- 



stable, and the disk is expected to exhibit the limit-cycle behavior (iKato et al.lll998l . Chaps. 
4 and 5). We have continued computations for several complete limit-cycles, and a repre- 
sentative cycle is illustrated in Figures [D - O which are for the time evolution of the radial 
distribution of the half-thickness of the disk H, temperature T, surface density S, effective 
optical depth t^jj = (2/3) (3x^/2 + + 1/tp), ratio of gas pressure to total pressure /3, 
and accretion rate m, respectively. Note that negative values of m signify an outflow in the 
radial direction, not in the vertical direction as the word 'outflow' in the literature usually 
means. 

The first panel of Figured] and the solid lines in Figures [2] - [6] show the disk just before 
the start of the cycle {t = Os). The disk is essentially in the SSD state, i.e., it is geometrically 
thin {H/r ^ 1) and optically thick {Tgff ^ 1) everywhere, its temperature has a peak at 
r ~ Gvg, its accretion rate is nearly constant in space, and its inner region (from ~ Sr^ 
to ~ 14rg) has (3 < 0.4 and is thermally unstable. Note that this configuration is not a 
stationary state and is with the diffusive viscosity prescription, so it is very different from 
the initial condition at the beginning of the computation, which is a stationary solution with 
the ap viscosity prescription. 

As the instability sets in (t = 2s, the second panel of Fig. [1] and the thin dashed lines 
in Figs. [2] -[6]), in the unstable region (r < 24rg) the temperature rises rapidly, the disk 
expands in the vertical direction, a very sharp spike appears in the surface density profile 
and accordingly in the optical depth and accretion rate profiles (exactly the stiff problem). 
The spikes move outwards with time, forming an expansion wave, heating the inner material 
and pushing it into the black hole, and perturbing the outer material to departure from 
the SSD state. The expansion region is in fact essentially in the state of slim disk, as it is 
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geometrically thick {H/r < 1), optically thick, very hot, and radiation pressure-supported 
{(3 < 0.4); and the front of the expansion wave forms the transition surface between the SSD 
state and the slim disk state. At t = 12s (the third panel of Fig. [T] and the thin dot-dashed 
lines in Figs. [2]-[n]), in the expansion region H and rh (negative, radial outflow) reach their 
maximum values, and the local rh (positive, inflow) exceeds 3 which is far above the initial 
value rh = 0.06 and is even well above the critical value m = 1. 

Once the wavefront has moved beyond the unstable region (r < 120rg), the expansion 
starts to weaken, the temperature drops in the innermost part of the disk and the material 
there deflates [t = 23s, the fourth panel of Fig. [T]and the thick dashed lines in Figs. [2] -[6]). 
Subsequently, deflation spreads out through the disk, and the disk consists of three different 
parts: the inner part is geometrically thin, with the temperature and surface density being 
lower than their values at t = Os; the middle part is what remains of the slim disk state; and 
the outer part is still basically in the SSD state (t = 27s, the fifth panel of Fig. [1] and the 
thick dot-dashed lines in Figs. [2] -[6]). 

The 'outburst' part of the cycle ends when it has proceeded on the thermal time-scale 
(t = 32s, the sixth panel of Fig. [T] and the dotted lines of Figs. [2]-[n]). What follows is a 
much slower process (on the viscous time-scale) of refilling and reheating of the inner part 
of the disk. Finally {t = 722s, the seventh panel of Fig. [1] and again the solid lines of Figs. 
[2] -[6]), the disk returns to essentially the same state as that at the beginning of the cycle. 
Then the thermal instability occurs again and a new cycle starts. 

The bolometric luminosity of the disk, obtained by integrating the radiated flux per unit 
area F~ over the disk at successive times, is drawn in Figure [7] for three complete cycles. 
The luminosity exhibits a burst with a duration of about 20 seconds and a quiescent phase 
lasting for the remaining about 700 seconds of the cycle. The amplitude of the variation 
is around two orders of magnitude, during the outburst a super-Eddington luminosity is 
realized. 

All these results obtained with our numerical method are similar to that of SMOl, 
not only in the sense that the limit-cycle behavior of thermally unstable accretion disks is 
confirmed, but also in the sense that the numerical solutions are of very good quality. In 
our computations we have been able to suppress all numerical instabihties and to remove all 
spurious oscillations, so that in our figures all curves are perfectly continuous and smooth 
and all spikes are well-resolved. 



What is new, however, is that we have also computed the Bernoulli functi on (i.e., the 
specific total energy) of the accreted matter that is expressed as (cf. Eq. [11.33] of lKato et al. 
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19981) 



B 



3(1-/3) + 



/3 



7 



(35) 



where 7 is the specific heat ratio and is taken to be 5/3. Figure [8] shows the quantity B 
obtained in the whole computational domain ranging from r ~ 2.5rg to r ~ lO^r^. It is clear 
that B has negative values for the whole spatial range (approaching to zero for very large r) 
and during the whole time of the cycle (in the figure the thick dot-dashed line for t = 27s 
and the dotted line for t = 32s are coincided with the solid line for t = Os). Note that in 
equation fl35l) the vertical kinetic energy 0.5V^ is included, and the gravitational energy is 
that for the surface of the disk. If the vertical kinetic energy is omitted and the gravitational 
energy is taken to be its equatorial plane value as in 1-dimensional models, then B will have 
ev en larger negative values. This result is in strong support of the analytical demonstration 
of lAbramowicz et al.l (120001 ) that accretion flows with not very strong viscosity {a < 0.1) 
have a negative Bernoulli function; and implies that outflows are unlikely to originate from 
thermally unstable accretion disks we consider here, because a positive B is a necessary, 
though not a sufficient, condition for the outflow production. 



4. Summary and Discussion 



We have introduced a numerical method for studies of thermally unstable accretion 
disks around black hole s, which is essent ially a combination of the standard one-domain 
pseudo-spectral met hod (IChan et al.ll2005l ) and the adaptive domain decomposition method 
( iBayliss et al.lll995l ). As a test of our method, for the case of moderate viscosity we have 
reproduced the best numerical results obtained previously by SMOl. Despite these simi- 
larities, we have made the following improvements over previous works in the numerical 
algorithm and concrete techniques, which have been proven to be effective in the practice of 
our computations. 

1. In applying the domain decomposition method to resolve the stiff problem, we de- 
velop a simple and useful connection technique to ensure a numerically stable continuity for 
the derivative of a physical quantity across the interface of tw o conjoint sub doma ins, i.e., 
equations ( fT3l) and dHj), instead of the connection technique in iBayliss et al.l ( 119951 ) that is 
seemingly complicated and was not explicitly explained. 



2. We construct a mapping function (eq. [H]) by adding the simple li near mapping 
functi on (eq. [16]) into the adaptive mapping function (eq. [17]) proposed by IBayliss et al. 



( 119951 ). so that the mapping between the Chebyshev-Gauss-Lobatto collocation points and 



the physical collocation points r^*'*, not only the mapping between two sets of collocation 
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points ffc and f^, is completed; and the adjustability of equation ffTTI) is kept to enable us 
to follow adaptively the region that is with the stiff problem and is shifting in space during 
time-evolution. 

3. For the time- integration, we use two complementary schemes, namely the third order 
TVD Runge-Kutta scheme and the third order backward-differentiation explicit scheme. The 
for mer scheme is pop ular in one-domain spectral methods and is essentially what was used 



by IChan et al.l (120051 ) . and the latter one can achieve the same accuracy and has advantage 



of lower CPU-time consumption. 



4. For the treatment o: 



boundary conditions, we notice that the spatial filter technique 



developed by I Chan et al.l (120051 ) for spectral methods is useful but is itself alone insuffi- 
cient, and the treatment traditionally used in finite difference methods is still needed to 
complement. We also find that once reasonable conditions are set at the outer boundary, 
our solutions behave themselves physically consistent close to the black hole horizon, and no 
inner boundary conditions are necessary as supplied by SMOl. 

5. We resolve the problem of spurious oscillations due to the absence of viscous stress 
tensor components in the basic equations in a way different from that of SMOl. SMOl 
introduced an artificial viscous term in the radial and vertical momentum equations. We 
instead improve the basic equations by including two viscous force terms Frr and in 
the radial momentum equation and a corresponding viscous heating term in the energy 
equati on, all these terms were alrea dy mentioned by the same authors of SMOl in an earlier 



paper (jSzuszkiewicz fc Miller! 119971 ). As for the vertical momentum equation, because of its 
crudeness in our 1.5-dimensional studies, we still adopt an artificial term whose explicit form 
is kindly provided by Szuszkiewicz & Miller and is unpublished. We obtain solutions at the 
same quality level as in SMOl, but we think that our treatment is probably more physical in 
some sense. In particular, any modification in the momentum equation ought to require a 
corresponding modification in the energy equation, otherwise the energy conservation would 
not be correctly described. 

Of these five improvements, we expect that the first two and the last one will be par- 
ticularly helpful for our subsequent studies of the strong viscosity case (a ~ 1). In this case 
the viscous heating becomes extremely huge, the 'outburst' of the disk due to the thermal 
instability is predicted to be more violent, and the Gibbs phenomenon related to the stiff 
problem can be even more serious than in the case of moderate viscosity studied in this 
paper. Our improved domain decomposition method is prepared to front these difficulties. 
As for another nettlesome problem that the absence of some viscous stress tensor compo- 
nents in 1- or 1.5-dimensional equations can also cause serious spurious oscillations, we think 
that, although in the moderate viscosity case they are equivalently effective as what were 
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made by SMOl, our modifications for both the radial momentum and energy equations will 
show their advantages in the strong viscosity case. In fa ct, the importance of the v iscous 



forces Frr and has long since been pointed out (e.g., iPapaloizou &: Stanley! Il986l ). We 
think that the inclusion of a heating term in the energy equation in accordance with these 
two forces will be not only consistent in physics, but also hopefully important in obtain- 
ing numerically stable solutions. With all these preparations made in this paper, we wish 
to achieve the goal to answer the question of the fate of thermally unstable black hole ac- 
cretion disks with very large values of a: d o these disks finally form sta ble and persistent 
SSD-I-ADAF configurations as suggested by lTakeuchi fc Mineshigd (Il998l ). or they also un- 
dergo limit-cycles, or something else? In view of the two facts that limit-cyclic luminosity 
variations, even though with seemingly very reliable theoretical warranties, are not usually 
observed for black hole systems (GRS 1915+105 remains the only one known); and that 
outfiows are already observed in many high energy astrophysical systems that are believed 
to be powered by black hole accretion, but are unlikely to originate from thermally unstable 
accretion disks we study here because of the negative Benoulli function of the matter in 
these disks, it will be definitely interesting if some behaviors other than the limit-cycle for 
non-stationary black hole accretion disks and/or the outfiow formation from these disks can 
be demonstrated theoretically. 
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A. Spectral Filtering and Boundary Conditions 



When applied to solve non-linear partial differential equations, a principle drawback of 
spectral methods is the aliasing error that can cause spurious oscillations at high frequen- 
cies. The spectral filtering is a special technique developed to filter out the high-f requency 
modes in each time-step to reduc e the aliasing error. As in IChan et al.l (120051 . see also 
Gottlieb fc Shul 119971 : |Peyretll2002l ). we use a exponential filter in spectral space as 



exp 



In. 



n 
N 



(Al) 



where e is the machine accuracy and 5 is a parameter that can be determined from the 
numerical practice. Then instead of u{rk), the filtered collocation values of the physical 
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quantity u{r) at a collocation point are given by 

N 

^NJ """" V 



crs { — ) u„, cos 



n=0 



(A2) 



As for the boundary conditions, the Dirichlet condition, u\boundary = Uq, or/and the 
Neumann condition, (du/dr) {boundary = Wg, are generally required. In order to avoid the 
appearance o f discontinuous ste p-functions at the boundaries that would cause the Gibbs 
Phenomenon, I Chan et al.l (l2005l ) introduced anot her filtering techn ique, namely the spatial 
filter. In our numerical algorithm, we either foUow lChan et al.l (l2005l ) to impose the boundary 
conditions by using the spatial filter, or directly change the value of a certain physical 
quantity to be its given value at the boundary point, depending on whether the boundary 
derivative of the quantity is or is not very small. The spatial filter is a monotonically 
decreasing filter 



h{r) = exp 



— I In el 



r — Tr 



(A3) 



for the outer boundary, and a monotonically increasing filter 



h{r) 



exp 



In el 



(A4) 



for the inner boundary. With this filter, the Dirichlet boundary condition can be imposed 
in each time-step as 

(A5) 



and the Neumann boundary condition can be imposed as 



du\' 
dr J ^ 



fduY , 



h{rk) + Mq, 



(A6) 



where the superscript n and subscript k denote the relevant values at the n-th time-level 
and the collocation point r^, respectively. 
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Fig. 2. — Evolution of the temperature of the matter in the disk during one full cycle. 
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3. — Evolution of the surface density of the disk. 
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Fig. 4.- 



Evolution of the effective optical deptfi of the disk. 
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5. — Evolution of the ratio of gas pressure to total pressure in the disk. 
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Fig. 6. — Evolution of the local accretion rate in the disk. 
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Fig. 7. — Variation of the bolometric luminosity of the disk during three full cycles. 
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Fig. 8. — Bernoulli function of the matter of the disk. Note that the horizonal scale is very 
different from that in Figs. [1] - [61 



